Systems of Linear Equations

In this notebook we leran various functions for solving and manipulating systems of linear equations.

Initialization

The main linear algebra routines are contained in scipy.linalg. Behind the scenes the calculations are done using BLAS and LAPACK libraries. Depending on how your version was compiled these can be highly optimized and even parallelized.

Notes:

1. Some functions also have implementations in NumPy. In general we are going to prefer the versions from scipy.linalg. In practice, both versions should give the same results except in extreme cases, such as dealing with (numerically) singular matrices.
2. We will represent both vectors and matrices using NumPy arrays. NumPy arrays are not vectors or matrices. A vector or matrix is a special object that has special properties. For example, we know what multiplying a matrix and a vector means, where as multiplying a two dimensional array and a one dimensional array is not uniquely defined. There is a matrix object in NumPy but we will not use those in this notebook. Due to the fact that we are using arrays instead of matrices this means we will need to use the function np.dot when we multiply a matrix and a vector. See below for details.

As always, begin by looking at the documentation. From scipy.linalg we see there are many routines. Here we will focus on some of those from "Basics" and those related to the LU decomposition.

Also as always we should look at the documentation for any and all functions that interest us.

Notice that there are a few functions related to the LU decomposition. We will learn why this is the case below.

Sample linear system

Consider a sample linear system of equations $$ \begin{align} 10 x_1 - 7 x_2 + x_3 &= 8, \\ -3 x_1 + 2 x_2 + 6 x_3 &= 4, \\ 5 x_1 - x_2 + 5 x_3 &= 6 . \end{align} $$ We can write this in the form $$\def\mat#1{\mathsf{#1}} A x = b, $$ where $$ A = \left( \begin{array}{rrr} 10 & -7 & 1 \\ -3 & 2 & 6 \\ 5 & -1 & 5 \end{array} \right) \quad\mathrm{and}\quad b = \left( \begin{array}{c} 8 \\ 4 \\ 6 \end{array} \right). $$

We can create these as arrays using np.array. Notice that we force the result to be a floating point array, not an integer array, by making any one of the entries a floating point number. NumPy tries to use the "simplest" data type when it creates an array. There are a few ways to force the type it chooses, this is one way.

Solving the System

To solve this system we can, not surprisingly, use solve.

Verifying the Solution

To verify the solution we can directly evaluate $A x$ and compare it to $b$. Of course this only verifies that the solve has worked correctly, it does not verify that we have entered $A$ and $b$ correctly. There are two steps, one is to actually perform the multiplication and the second is to compare the $A x$ to $b$.

Matrix Multiplication

To multiply matrices we do not use '*' . The usual multiplication is component-wise multiplication. This is not the same as matrix multiplication.

To see this, consider using the usual multiplication:

Notice what happened: the first component of $b$ multiplied the first column of $A$, the second component of $b$ multiplied the second column of $A$, and similarly for the third. This is what we mean by component-wise multiplication.

If we want to perform matrix multiplication we need to use a function, np.dot. Consider

Notice that this returns a vector and that vector should be $b$.

Comparing calculations

This is a small system so we can easily see that all values are nearly zero. We do not want to check that $A x$ is exactly equal to $b$ (in general this does not work numerically), however, we can check that the values are "close enough". There are a few ways to do this, here we use np.allclose. (At this point you may be checking the documentation! Here we are using the default settings, but could adjust the tolerances for a more stringent test.)

Factoring the system

The work required to solve the system of equations mostly involves manipulating the matrix, $A$, and then performing the same manipulations on the right hand side of the equations, $b$. We could instead have many right hand sides (a two dimensional array with multiple columns, one for each set of values for which we want to find a solution). This is handled by solve. Alternatively, we can decompose $A$ into pieces that encode the required manipulations using the LU decomposition. The decomposition only needs to be performed once, it can then be applied whenever needed. Finally, for numerical stability we should also pivot the matrix with permutation matrix $P$. LU decomposition with pivoting is provided by scipy.linalg.lu.

By default we are given the permutation matrix $P$ and the lower and upper triangular matrices that, when recombined, produce $\mat A$. This function is good when we want to see the decomposition in a form easy for us to read. If we want to use the decomposition for solving systems of linear equations the information could be stored in a more efficient form for the computer's use.

Note that our choice of using generic arrays for storing matrices comes at a cost when we want to multiply many of them together. Here we have had to use nested calls to np.dot. This can be slightly streamlined, but remains tedious. It is just a price we must pay for our choice. In practice it is usually a rather small price.

Getting back to solving a system of equations, we can use scipy.linalg.lu_factor to decompose in a form more useful for the computer.

Here $L$ and $U$ have been merged into a single matrix (notice in the form above that $L$ had ones along the diagonal so they do not need to be stored) and the permutations are represented by a one dimensional array instead a matrix. This is much more memory efficient, but is also much harder for us to read. Even so, this can be used in scipy.linalg.lu_solve. In fact, the tuple returned here is exactly what needs to be provided as the first argument to that function.

If $A=LU$, then $Ax=LUx=b$. Thus, first solve $Ly=b$, for $y$, then solve $Ux=y$ for $x$. Both solutions are found by back-substitution.

With permutation matrix $P$, such that $PA=LU$ the decomposition is sometimes called LUP decomposition, so $LUx=Pb$. In this case the solution is also done in two logical steps: (1) solve the equation $Ly=Pb$ for $y$ and (2) solve the equation $Ux=y$ for $x$.

Slightly Larger System

As an example of a slightly larger system and one where we want to find solutions with different right hand sides consider $$\mat A = \left( \begin{array}{rrrr} 2 & -3 & 1 & 3 \\ 1 & 4 & -3 & -3 \\ 5 & 3 & -1 & -1 \\ 3 & -6 & -3 & 1 \end{array} \right), \quad \vec b_1 = \left( \begin{array}{r} -4 \\ 1 \\ 8 \\ -5 \end{array} \right), \quad \vec b_2 = \left( \begin{array}{r} -10 \\ 0 \\ -3 \\ -24 \end{array} \right). $$

We can solve this directly using solve.

Alternatively we can use an LU decomposition. Notice the decomposition is only performed once.

Finally, we could have solved for both right hand sides at once. To do this we need to combine b1 and b2. We do this below.

If we try to use this with solve it fails!

We get an incompatible dimensions error because b should have the same number of rows as A, not the same number of columns. This is easy to fix, just take the transpose of b. There is shorthand for doing this as shown below.

With this we find:

Each column of the result is the solution for each column in the array b.

Note: This notebook is adapted from http://www.phys.cwru.edu/courses/p250/